Path: hard-problems-formulations / minimum-set-cover.ipynb

Views: ¹⁰⁰
Image: ubuntu2004

Kernel: .venv

MINIMUM SET COVER / Минимальное покрытие множеств

In [1]:

import random
from IPython.core.display import SVG
import pyomo.environ as pyo
from pysat.solvers import Solver
from pysat.formula import CNF 
import py_svg_combinatorics as psc
from ipywidgets import widgets, HBox
from collections import Counter
from pprint import pprint
from random import randint
import numpy as np
from IPython.display import IFrame
import IPython
from copy import copy
import os
from pathlib import Path
nbname = ''
try:
    nbname = __vsc_ipynb_file__
except:
    if 'COCALC_JUPYTER_FILENAME' in os.environ:
        nbname = os.environ['COCALC_JUPYTER_FILENAME']    
title_ = Path(nbname).stem.replace('-', '_').title()    
IFrame(f'https://discopal.ispras.ru/index.php?title=Hardprob/{title_}&useskin=cleanmonobook', width=1280, height=300)

In [2]:

import random
from IPython.core.display import SVG
import pyomo.environ as pyo
from pysat.solvers import Solver
from pysat.formula import CNF 
import py_svg_combinatorics as psc
from ipywidgets import widgets, HBox
from collections import Counter
from pprint import pprint
import numpy as np

Постановка задачи

Представим, что мы имеем конечное множество «U» (universum) и «S» — список подмножеств множества «U».

'''Покрытием''' называют семейство $\mathcal{C}\subseteq\mathcal{S}$ наименьшей мощности, объединением которых является $\mathcal{U}$ . В случае постановки вопроса '''о разрешении''' на вход подаётся пара $(\mathcal{U},\mathcal{S})$ и целое число k; вопросом является существование покрывающего множества мощности k (или менее).

Оптимизационная версия задачи, '''минимальное покрытие множеств''', задаёт вопрос о минимальном числе множеств из «S» покрывающих «U».

https://en.wikipedia.org/wiki/Set_cover_problem

In [3]:

sets = [
        [f"{i:02}" for i in range(1,11)],
#        [f"{i:02}" for i in range(11,21)],
        [f"{2*i:02}" for i in range(1,10)]+['01'],
        [f"{i:02}" for i in np.random.choice(np.arange(1,20), 1)],
        [f"{i:02}" for i in np.random.choice(np.arange(1,20), 2)],
        [f"{i:02}" for i in np.random.choice(np.arange(1,20), 3)],
        [f"{i:02}" for i in np.random.choice(np.arange(1,20), 4)],
        [f"{i:02}" for i in np.random.choice(np.arange(1,20), 5)],
        [f"{i:02}" for i in np.random.choice(np.arange(1,20), 6)],
        [f"{i:02}" for i in np.random.choice(np.arange(1,20), 7)],
        [f"{i:02}" for i in np.random.choice(np.arange(1,20), 8)],
        [f"{i:02}" for i in np.random.choice(np.arange(1,20), 9)],
]

In [4]:

svg = psc.subsets2svg(sets)
SVG(data=svg)

Реализация в Pyomo

In [5]:

def print_solution(m):
    for v in m.component_data_objects(pyo.Var):
        if v.value and v.value > 0:
            print(str(v), v.value)

In [6]:

def get_model(sets):
    sets_ = sets
    if isinstance(sets_, dict):
        sets_ = [val for _, val in sets_.items()]

    model = pyo.ConcreteModel()
    model.S = [set(str(item) for item in sublist) for sublist in sets_] 
    model.m = len(model.S)
    model.J = range(model.m)
    model.U = sorted(set([str(item) for sublist in sets_ for item in sublist]))
    model.n = len(model.U)

    model.x = pyo.Var(model.J, domain=pyo.Binary)
    model.set_count = pyo.Objective(expr = sum( model.x[j] for j in model.J), sense=pyo.minimize)

    @model.Constraint(model.U)
    def каждый_элемент_в_одном_множестве(m, u):
        return sum(m.x[j] for j in m.J if u in m.S[j]) >= 1

    return model


m = get_model(sets)
m.U, m.S

(['01',
  '02',
  '03',
  '04',
  '05',
  '06',
  '07',
  '08',
  '09',
  '10',
  '11',
  '12',
  '13',
  '14',
  '15',
  '16',
  '17',
  '18'],
 [{'01', '02', '03', '04', '05', '06', '07', '08', '09', '10'},
  {'01', '02', '04', '06', '08', '10', '12', '14', '16', '18'},
  {'03'},
  {'04', '11'},
  {'09', '11', '16'},
  {'03', '04', '17'},
  {'02', '03', '07', '10', '18'},
  {'09', '10', '12', '13', '15', '16'},
  {'03', '05', '06', '08', '12', '17'},
  {'01', '04', '06', '07', '11', '14', '18'},
  {'01', '04', '05', '12', '13', '15', '17'}])

In [7]:

m.pprint()

2 Set Declarations
    x_index : Size=1, Index=None, Ordered=Insertion
        Key  : Dimen : Domain : Size : Members
        None :     1 :    Any :   11 : {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10}
    каждый_элемент_в_одном_множестве_index : Size=1, Index=None, Ordered=Insertion
        Key  : Dimen : Domain : Size : Members
        None :     1 :    Any :   18 : {'01', '02', '03', '04', '05', '06', '07', '08', '09', '10', '11', '12', '13', '14', '15', '16', '17', '18'}

1 Var Declarations
    x : Size=11, Index=x_index
        Key : Lower : Value : Upper : Fixed : Stale : Domain
          0 :     0 :  None :     1 : False :  True : Binary
          1 :     0 :  None :     1 : False :  True : Binary
          2 :     0 :  None :     1 : False :  True : Binary
          3 :     0 :  None :     1 : False :  True : Binary
          4 :     0 :  None :     1 : False :  True : Binary
          5 :     0 :  None :     1 : False :  True : Binary
          6 :     0 :  None :     1 : False :  True : Binary
          7 :     0 :  None :     1 : False :  True : Binary
          8 :     0 :  None :     1 : False :  True : Binary
          9 :     0 :  None :     1 : False :  True : Binary
         10 :     0 :  None :     1 : False :  True : Binary

1 Objective Declarations
    set_count : Size=1, Index=None, Active=True
        Key  : Active : Sense    : Expression
        None :   True : minimize : x[0] + x[1] + x[2] + x[3] + x[4] + x[5] + x[6] + x[7] + x[8] + x[9] + x[10]

1 Constraint Declarations
    каждый_элемент_в_одном_множестве : Size=18, Index=каждый_элемент_в_одном_множестве_index, Active=True
        Key : Lower : Body                                     : Upper : Active
         01 :   1.0 :               x[0] + x[1] + x[9] + x[10] :  +Inf :   True
         02 :   1.0 :                       x[0] + x[1] + x[6] :  +Inf :   True
         03 :   1.0 :         x[0] + x[2] + x[5] + x[6] + x[8] :  +Inf :   True
         04 :   1.0 : x[0] + x[1] + x[3] + x[5] + x[9] + x[10] :  +Inf :   True
         05 :   1.0 :                      x[0] + x[8] + x[10] :  +Inf :   True
         06 :   1.0 :                x[0] + x[1] + x[8] + x[9] :  +Inf :   True
         07 :   1.0 :                       x[0] + x[6] + x[9] :  +Inf :   True
         08 :   1.0 :                       x[0] + x[1] + x[8] :  +Inf :   True
         09 :   1.0 :                       x[0] + x[4] + x[7] :  +Inf :   True
         10 :   1.0 :                x[0] + x[1] + x[6] + x[7] :  +Inf :   True
         11 :   1.0 :                       x[3] + x[4] + x[9] :  +Inf :   True
         12 :   1.0 :               x[1] + x[7] + x[8] + x[10] :  +Inf :   True
         13 :   1.0 :                             x[7] + x[10] :  +Inf :   True
         14 :   1.0 :                              x[1] + x[9] :  +Inf :   True
         15 :   1.0 :                             x[7] + x[10] :  +Inf :   True
         16 :   1.0 :                       x[1] + x[4] + x[7] :  +Inf :   True
         17 :   1.0 :                      x[5] + x[8] + x[10] :  +Inf :   True
         18 :   1.0 :                       x[1] + x[6] + x[9] :  +Inf :   True

5 Declarations: x_index x set_count каждый_элемент_в_одном_множестве_index каждый_элемент_в_одном_множестве

In [8]:


ilp_solver = pyo.SolverFactory('cbc')
ilp_solver.solve(m).write()
print_solution(m)

# ==========================================================
# = Solver Results                                         =
# ==========================================================
# ----------------------------------------------------------
#   Problem Information
# ----------------------------------------------------------
Problem: 
- Name: unknown
  Lower bound: 4.0
  Upper bound: 4.0
  Number of objectives: 1
  Number of constraints: 17
  Number of variables: 11
  Number of binary variables: 11
  Number of integer variables: 11
  Number of nonzeros: 11
  Sense: minimize
# ----------------------------------------------------------
#   Solver Information
# ----------------------------------------------------------
Solver: 
- Status: ok
  User time: -1.0
  System time: 0.01
  Wallclock time: 0.01
  Termination condition: optimal
  Termination message: Model was solved to optimality (subject to tolerances), and an optimal solution is available.
  Statistics: 
    Branch and bound: 
      Number of bounded subproblems: 0
      Number of created subproblems: 0
    Black box: 
      Number of iterations: 0
  Error rc: 0
  Time: 0.5427391529083252
# ----------------------------------------------------------
#   Solution Information
# ----------------------------------------------------------
Solution: 
- number of solutions: 0
  number of solutions displayed: 0
x[0] 1.0
x[7] 1.0
x[9] 1.0
x[10] 1.0

In [9]:

print_solution(m)

x[0] 1.0
x[7] 1.0
x[9] 1.0
x[10] 1.0

In [10]:

selected = [j for j in m.J if m.x[j].value > 0]
selected

[0, 7, 9, 10]

In [11]:

SVG(psc.subsets2svg(sets, selected))

In [ ]:

In [12]:

cnf3 = CNF(from_clauses=psc.rand3cnf(5, 3))
cnf3.clauses

[(-3, -1, -2), (2, 3, -1), (-2, 1, 3), (1, 2, 3), (3, -2, -1)]

In [13]:

from pysat.solvers import Solver

solver = Solver(bootstrap_with=cnf3)
res = solver.solve()
res

True

In [14]:

print(solver.get_model())

[1, -2, 3]

| | |

In [15]:

def reduce_3sat_to_msc(phi):
    # phi: a 3SAT formula with n variables and m clauses

    # create the universe U
    U = set()
    for clause in phi:
        for literal in clause:
            U.add(abs(literal))

    # create the collection of subsets S
    S = []
    for clause in phi:
        Sj = set()
        for literal in clause:
            v = abs(literal)
            if literal > 0:
                Sj.add(v * 2 - 1)  # add v_1
            else:
                Sj.add(v * 2)  # add v_0
        S.append(Sj)

    # set k = m
    k = len(phi)

    return U, S, k

In [16]:

U, S, k = reduce_3sat_to_msc(cnf3.clauses)
U, S, k

({1, 2, 3}, [{2, 4, 6}, {2, 3, 5}, {1, 4, 5}, {1, 3, 5}, {2, 4, 5}], 5)

In [17]:

cnf3.clauses

[(-3, -1, -2), (2, 3, -1), (-2, 1, 3), (1, 2, 3), (3, -2, -1)]

In [18]:

SVG(psc.subsets2svg(S))

In [19]:

print(list(solver.enum_models()))

[[1, -2, 3], [-1, -2, 3], [-1, 2, 3]]

In [20]:

from collections import defaultdict
def cnf32setcovering(cnf):
    '''
    Будем считать, что на входе 3SAT формула без повторов литералов в дизьюнкциях и т.п.
    '''
    n = cnf.nv
    m = len(cnf.clauses)
    S = {}

    for j in range(1, n+1):
        S[f'x_{{{j}}}=0'] = [f'x_{{{j}}}']
        S[f'x_{{{j}}}=1'] = [f'x_{{{j}}}']

    for i in range(m):
        U_clause = f'''C_{{{i}}}'''
        clause = cnf.clauses[i]
        for lit in clause:
            j = abs(lit)
            needvar = f'x_{{{j}}}=1'
            if lit < 0:
                needvar = f'x_{{{j}}}=0'
            S[needvar].append(U_clause)

    return S

In [21]:

cnf3 = CNF(from_clauses=psc.rand3cnf(1, 3))
cnf3.clauses

[(3, -1, -2)]

In [22]:

sets = cnf32setcovering(cnf3)
pprint(sets)

{'x_{1}=0': ['x_{1}', 'C_{0}'],
 'x_{1}=1': ['x_{1}'],
 'x_{2}=0': ['x_{2}', 'C_{0}'],
 'x_{2}=1': ['x_{2}'],
 'x_{3}=0': ['x_{3}'],
 'x_{3}=1': ['x_{3}', 'C_{0}']}

In [23]:

SVG(psc.subsets2svg(sets))

In [24]:

m = get_model(sets)
m.U, m.S

(['C_{0}', 'x_{1}', 'x_{2}', 'x_{3}'],
 [{'C_{0}', 'x_{1}'},
  {'x_{1}'},
  {'C_{0}', 'x_{2}'},
  {'x_{2}'},
  {'x_{3}'},
  {'C_{0}', 'x_{3}'}])

In [25]:

# ilp_solver = pyo.SolverFactory('cbc')
ilp_solver.solve(m).write()
print_solution(m)

# ==========================================================
# = Solver Results                                         =
# ==========================================================
# ----------------------------------------------------------
#   Problem Information
# ----------------------------------------------------------
Problem: 
- Name: unknown
  Lower bound: 3.0
  Upper bound: 3.0
  Number of objectives: 1
  Number of constraints: 4
  Number of variables: 6
  Number of binary variables: 6
  Number of integer variables: 6
  Number of nonzeros: 6
  Sense: minimize
# ----------------------------------------------------------
#   Solver Information
# ----------------------------------------------------------
Solver: 
- Status: ok
  User time: -1.0
  System time: 0.01
  Wallclock time: 0.01
  Termination condition: optimal
  Termination message: Model was solved to optimality (subject to tolerances), and an optimal solution is available.
  Statistics: 
    Branch and bound: 
      Number of bounded subproblems: 0
      Number of created subproblems: 0
    Black box: 
      Number of iterations: 0
  Error rc: 0
  Time: 0.07667303085327148
# ----------------------------------------------------------
#   Solution Information
# ----------------------------------------------------------
Solution: 
- number of solutions: 0
  number of solutions displayed: 0
x[1] 1.0
x[3] 1.0
x[5] 1.0

In [26]:

selected = [j for j in m.J if m.x[j].value > 0]
selected

[1, 3, 5]

In [27]:

cnf3.clauses

[(3, -1, -2)]

In [28]:

pprint(sets)

{'x_{1}=0': ['x_{1}', 'C_{0}'],
 'x_{1}=1': ['x_{1}'],
 'x_{2}=0': ['x_{2}', 'C_{0}'],
 'x_{2}=1': ['x_{2}'],
 'x_{3}=0': ['x_{3}'],
 'x_{3}=1': ['x_{3}', 'C_{0}']}

In [29]:

SVG(psc.subsets2svg(sets, selected))

In [ ]:

In [30]:

cnf3 = CNF(from_clauses=psc.rand3cnf(2, 3))
cnf3.clauses

[(-3, 2, 1), (3, 1, 2)]

In [31]:

sets = cnf32setcovering(cnf3)
pprint(sets)

{'x_{1}=0': ['x_{1}'],
 'x_{1}=1': ['x_{1}', 'C_{0}', 'C_{1}'],
 'x_{2}=0': ['x_{2}'],
 'x_{2}=1': ['x_{2}', 'C_{0}', 'C_{1}'],
 'x_{3}=0': ['x_{3}', 'C_{0}'],
 'x_{3}=1': ['x_{3}', 'C_{1}']}

In [32]:

SVG(psc.subsets2svg(sets))

In [33]:

m = get_model(sets)
m.U, m.S

(['C_{0}', 'C_{1}', 'x_{1}', 'x_{2}', 'x_{3}'],
 [{'x_{1}'},
  {'C_{0}', 'C_{1}', 'x_{1}'},
  {'x_{2}'},
  {'C_{0}', 'C_{1}', 'x_{2}'},
  {'C_{0}', 'x_{3}'},
  {'C_{1}', 'x_{3}'}])

In [34]:

ilp_solver.solve(m).write()
print_solution(m)

# ==========================================================
# = Solver Results                                         =
# ==========================================================
# ----------------------------------------------------------
#   Problem Information
# ----------------------------------------------------------
Problem: 
- Name: unknown
  Lower bound: 3.0
  Upper bound: 3.0
  Number of objectives: 1
  Number of constraints: 5
  Number of variables: 6
  Number of binary variables: 6
  Number of integer variables: 6
  Number of nonzeros: 6
  Sense: minimize
# ----------------------------------------------------------
#   Solver Information
# ----------------------------------------------------------
Solver: 
- Status: ok
  User time: -1.0
  System time: 0.01
  Wallclock time: 0.0
  Termination condition: optimal
  Termination message: Model was solved to optimality (subject to tolerances), and an optimal solution is available.
  Statistics: 
    Branch and bound: 
      Number of bounded subproblems: 0
      Number of created subproblems: 0
    Black box: 
      Number of iterations: 0
  Error rc: 0
  Time: 0.040219783782958984
# ----------------------------------------------------------
#   Solution Information
# ----------------------------------------------------------
Solution: 
- number of solutions: 0
  number of solutions displayed: 0
x[1] 1.0
x[2] 1.0
x[5] 1.0

In [35]:

selected = [j for j in m.J if m.x[j].value > 0]
selected

[1, 2, 5]

In [36]:

cnf3.clauses

[(-3, 2, 1), (3, 1, 2)]

In [37]:

pprint(sets)

{'x_{1}=0': ['x_{1}'],
 'x_{1}=1': ['x_{1}', 'C_{0}', 'C_{1}'],
 'x_{2}=0': ['x_{2}'],
 'x_{2}=1': ['x_{2}', 'C_{0}', 'C_{1}'],
 'x_{3}=0': ['x_{3}', 'C_{0}'],
 'x_{3}=1': ['x_{3}', 'C_{1}']}

In [38]:

SVG(psc.subsets2svg(sets, selected))

In [39]:

m.pprint()
m.U

2 Set Declarations
    x_index : Size=1, Index=None, Ordered=Insertion
        Key  : Dimen : Domain : Size : Members
        None :     1 :    Any :    6 : {0, 1, 2, 3, 4, 5}
    каждый_элемент_в_одном_множестве_index : Size=1, Index=None, Ordered=Insertion
        Key  : Dimen : Domain : Size : Members
        None :     1 :    Any :    5 : {'C_{0}', 'C_{1}', 'x_{1}', 'x_{2}', 'x_{3}'}

1 Var Declarations
    x : Size=6, Index=x_index
        Key : Lower : Value : Upper : Fixed : Stale : Domain
          0 :     0 :   0.0 :     1 : False : False : Binary
          1 :     0 :   1.0 :     1 : False : False : Binary
          2 :     0 :   1.0 :     1 : False : False : Binary
          3 :     0 :   0.0 :     1 : False : False : Binary
          4 :     0 :   0.0 :     1 : False : False : Binary
          5 :     0 :   1.0 :     1 : False : False : Binary

1 Objective Declarations
    set_count : Size=1, Index=None, Active=True
        Key  : Active : Sense    : Expression
        None :   True : minimize : x[0] + x[1] + x[2] + x[3] + x[4] + x[5]

1 Constraint Declarations
    каждый_элемент_в_одном_множестве : Size=5, Index=каждый_элемент_в_одном_множестве_index, Active=True
        Key   : Lower : Body               : Upper : Active
        C_{0} :   1.0 : x[1] + x[3] + x[4] :  +Inf :   True
        C_{1} :   1.0 : x[1] + x[3] + x[5] :  +Inf :   True
        x_{1} :   1.0 :        x[0] + x[1] :  +Inf :   True
        x_{2} :   1.0 :        x[2] + x[3] :  +Inf :   True
        x_{3} :   1.0 :        x[4] + x[5] :  +Inf :   True

5 Declarations: x_index x set_count каждый_элемент_в_одном_множестве_index каждый_элемент_в_одном_множестве

['C_{0}', 'C_{1}', 'x_{1}', 'x_{2}', 'x_{3}']

In [ ]: